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^ | Abstract. The entangled quantum state of a photon pair propagating through 

atmospheric turbulence suffers decay of entanglement due to the scintillation it 
experiences. Here we investigate the robustness against this decay for different qutrit 
states. We use an inhnitesimal-propagation equation to obtain the density matrix as a 
function of the propagation distance and we use the tangle to quantify the entanglement 
between a pair of qutrits. The evolution of various initial states as they propagate 
through turbulence is considered. Using optimization of the initial parameters, we 
. obtain expressions for bipartite qutrit states that retain their initial entanglement 

longer than the initially maximally entangled states. 
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1. Introduction 

Free-space quantum communication is one of the major components in the new quantum 
information technology revolution. A challenge that confronts free-space quantum 
communication is the decay of entanglement that entangled photons experience due 
to scintillation while propagating through a turbulent atmosphere. The approaches to 
overcome this challenge are either based on methods to correct the output optical field, 
such as using adaptive optics [T|, or on choosing an input optical field that is to some 
extent robust against entanglement decay in atmospheric turbulence. Here we consider 
the latter approach. In other words, we investigate to what extent a quantum state can 
be optimized so that it will retain as much of its initial entanglement as possible, while 
propagating through turbulence. We'll refer to such an optimized quantum state as a 
robust state. 

In this paper we consider bipartite qutrit states in the Laguerre-Gaussian (LG) 
modal basis, which we restrict to the three elements that have a radial modal index of 
zero (p = 0) and azimuthal modal indices of i = 1, 0, —1, which gives a three-dimensional 
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Hilbert space per photon H 3 . The orbital angular momentum (OAM) associated with 
an LG mode is proportional to I. If, instead of choosing \t\ < 1, we use larger values of 
I (higher OAM), the entanglement is expected to last longer [2j. This follows from the 
observation that the coupling from one OAM mode into another is stronger the closer 
the two OAM values are to each other [31 H]. 

Here the infinitesimal-propagation equation (IPE) [5j is used to calculate the z- 
dependence of the corresponding density matrix. We assume that the two photons 
propagate through different uncorrelated regions of atmospheric turbulence. The 
turbulence is modeled by the Kolmogorov power spectral density [6]. 

We quantify the entanglement between the two photons by the tangle r IE] as a 
function of z. The tangle is equal to the square of the concurrence for pure states and 
gives a lower bound for the square of the concurrence of mixed quantum states. To find 
the most robust qutrit states, we optimize the tangle at a given propagation distance 
z > 0. The result is expressed in terms of the parameters for the initial pure state that 
will give the maximum tangle at that propagation distance. 

It has been shown that if only one partie of a bipartite (or multipartite) state passes 
through a dissipative channel, the final amount of entanglement is proportional to the 
initial entanglement of the state [9j [TOj [HI [12] . For such cases the initially maximally 
entangled states are also the most robust states. This argument has been extended 
to the case where all parties of an entangled state pass through dissipative channels, 
provided that these channels are trace preserving [13]. However, if both parties of an 
entangled bipartite state pass through uncorrelated dissipative channels that are not 
trace preserving, only an upper bound exists for the evolution of the entanglement 
[HJ [TOl EH]. Here we'll show in particular that for two photons propagating through a 
turbulent atmosphere, corresponding to two qutrit states passing through uncorrelated 
dissipative channels that are not trace preserving, the initially maximally entangled 
state is not also the most robust state. 

2. Initial state 



The most general pure bipartite qutrit state, defined in terms of our Hilbert space 
U 3 ® U 3 is 



where c mj „ represents complex coefficients and \m) and |n) are the OAM eigenstates for 
the respective subsystems, with m,n 6 {1,0, —1}. The normalization of the initial state 
implies that 



We use a parameterization that implicitly obeys this normalization condition. It also 
removes an overall phase factor and incorporates the symmetry associated with an 




(1) 



m,n 




(2) 



m,n 



Robust states 



3 



interchange of the I = 1 and £ = — 1 states. This parameterization is given by 
ci,i = sin(£; a ) sin(£; b ) sin(A; d ) sin(/c h ) exp[i(g a + g b + q d + g d )] 
ci _i = cos(£; a ) sin(fc b ) sin(/c d ) sin(fc h ) exp[i(-g a + g b + g d + g h )] 
c_ X) i = sin(& c ) cos(& b ) sin(A; d ) sin(fc h ) exp[i(g c - q h + q d + q h )} 
C-1,-1 = cos(A; c ) cos(fcb) sin(fc d ) sm(k h ) exp[i(-g c - q h + q d + q h )] 
c ,o = cos(/c h ) exp(-ig h ) (3) 
ci j0 = sin(/c e ) sin(fcf) cos(/c d ) sin(fc h ) exp[i(g e + g f - g d + gh)] 
c_i,q = cos(A; e ) sin(/c f ) cos(fc d ) sin(fc h ) exp[i(-g c + g f - g d + g h )] 
c 0) i = sin(fcg) cos(/c f ) cos(/c d ) sin(fc h ) exp[i(g g - q f — q d + g h )] 
Co _i = cos(A; g ) cos(/c f ) cos(A; d ) sin(fc h ) exp[i(-g g - g f — g d + gh)]- 

Due to the normalization condition and the removal of the overall phase we have 16 real 
parameters: 8 angles (fc a , fch) an d 8 phases (g a , g^). 

The density matrix p(z) is obtained as a function of the propagation distance by 
solving the IPE [5] for an initial density matrix defined as p(0) = where is 

given in (pQ), using the coefficients defined in ()3]). Although, we assume that the initial 
state (OQ) is a pure state, the scintillation process causes it to become mixed during 
propagation, which necessitates a density matrix approach. 



3. The IPE 



The IPE [5] is a set of coupled first order differential equations that describes the 
evolution of the density matrix for a biphoton in the OAM basis as a function of the 
propagation distance z through a turbulent medium. It represents a multiple phase 
screen approach, as opposed to the single phase screen approach [3j. As a result it can 
simulate both phase and intensity fluctuations, whereas the single phase screen approach 
can only simulate phase fluctuations. 

In the IPE the turbulence model is specified in terms of a power spectral density. 
Here we'll use the Kolmogorov power spectral density [6], given by 

$(£;) = 0.033 Clk- n/3 , (4) 

where C\ is the refractive index structure constant, which quantifies the strength of the 
turbulence. In the solution of the density matrix, obtained from the IPE the strength 
of the turbulence is contained in a dimensionless parameter defined by 

^3/2(^2 11/3 

= 71 L 'n w /r\ 

6r(2/3)A 3 ' { ' 

where wq is the radius of the Gaussian envelope at the waist of the beam, A is the 
wavelength of the photons and is the Rayleigh range, defined by 

ZR = ^. (6) 
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For the three-dimensional bipartite case the density matrix is a 9 x 9 matrix. 
Therefore, the IPE for this density matrix is given as 81 first order differential equations 
that are divided into several sets of coupled equations and a few uncoupled equations. To 
solve the sets of coupled differential equations we use a perturbative approach, exploiting 
the fact that the terms that couple the different equations are multiplied by constants 
that are invariably an order of magnitude smaller that the uncoupled terms. 

After solving the differential equations to obtain the density matrix as a function 
of the propagation distance, one finds that its dependence on the propagation distance 
is governed by two functions that are defined by the following integrals 

Z(t) = Al + r 2 ) 5 / 6 dr (7) 



and 



H(t) 
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where we define a normalized propagation distance t = z/z-r. For convenience we 
separate H{t) into its real and imaginary parts H(t) = H T {t) + \H{(t). These integrals 
can be solved to give expressions in terms of hypergeometric functions 
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where 2Fi([ - ]> [•],*) represents Barnes's extended hypergeometric functions 
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Figure 1: Full and approximate solutions of the integrals for Z(t), H r (t) and H\(t) as a 
function of the normalized propagation distance t. 

In figure [1] we show the exact solutions for Z(t), H r (t) and H\(i) in comparison 
with their approximations, indicating that if propagation is limited to t < 1/3, one 
can use the approximate solutions. The requirement that t < 1/3 implies that the 
entanglement decays quickly, which calls for strong turbulence conditions. In what 
follows we'll consider both the cases where the entanglement decays at t < 1/3 (strong 
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turbulence), using the approximate expressions in (fTD|) and (TTTT) . and the cases where 
the entanglement only decays at t > 1/3 (weak turbulence), using the full expressions 
in (ED, flE]) and AH]). 

4. Trace 

The solution that is obtained from the IPE represents a truncated density matrix, 
because only three of the infinite number of basis elements are retained. This is 
experimentally analogous to post selection [E] of measurement results in the three 
dimensional subspace. During propagation through turbulence the energy in the initial 
modes is scattered into higher order modes that are not contained in the truncated 
density matrix. As a result the density matrix is subnormalized — its trace is smaller 
that 1. This, in turn, causes the trace of the truncated density matrix to decrease as a 
function of the propagation distance. Moreover, the function of the trace depends on the 
initial state. However, due to the symmetries among the chosen OAM basis elements, 
the trace function only contains the parameters k^ and k d - 




Figure 2: Trace evolution of the truncated density matrix for strong turbulence (a) 
and for weak turbulence (b). Each graph contains three curves: solid (blue) line with 
fch = 0, short dashed (red) line with = ir/2, k d = and long dashed (green) line with 
kh = kd = 7r/2. 

The evolution of the trace as a function of t is plotted in figure |2] for three extreme 
cases: k h = 0, which is independent of k d ; k h = ir/2, k d = 0; and k h = k d = n/2. These 
cases are considered for strong turbulence (cr = 25) in figure l2al and for weak turbulence 
(a = 0.25) in figure |2bl The plots show that states with I = scatter slower than 
those with i = ±1. In fact, the slowest decay of the trace is obtained for k^ = 0, which 
represents a state consisting of only 1 , 0) and, therefore, is completely separable. 
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5. Tangle 

For pure two-dimensional bipartite systems the concurrence [16] is a suitable measure of 
entanglement, but for mixed states in higher dimensions it is computationally demanding 
to calculate the concurrence directly through the construction of a convex roof [17J. 
Instead, we use the tangle r{p}, which is equal to the square of the concurrence for 
pure states and gives a lower bound for that of mixed quantum states [7]. The tangle 
is calculated via the purities of the (reduced) density matrices 

r{p} = 2 tr{p 2 } - tr{pi} - tr{p 2 B }, (12) 

where pa and p-Q are the respective reduced density matrices of the two subsystems. For 
a maximally entangled state r = r max = 2(d — l)/d, where d is the dimension of the 
subsystems' Hilbert spaces and for a separable state r = r min = 0. 



5.1. Bell states 

In the two-dimensional case, where diametric azimuthal indices I = ±1 are used, it can 
be shown that the entanglement of the four maximally entangled Bell states, quantified 
by the concurrence, all decay equally [5]. This can be understood as a result of the 
fact that the scattering is symmetric with respect to I = 1 and I = — 1. In the three- 
dimensional case there is a third possible azimuthal index, which does not share this 
symmetry. In our case we chose the third index to be i = 0, and we quantify the 
entanglement by the tangle. As a result, we find that the 12 Bell states in the two- 
dimensional subspaces of our three-dimensional Hilbert space form three sets based on 
their decay curves 

Set 1 = {|<&+) , !<&+_!> , 1^) , 1$^)} (13) 
Set 2 = {jtf+A , , , (14) 

Set 3 = , |$ w > , !*+_!) , (15) 

where 

|*±) = -^(|r,r>±M) (16) 

|*±> = -^(|r,a>±| aj r» (17) 

are the Bell states of the two-dimensional subspaces spanned by r and s. The three sets 
of subspace Bell states give three different decay curves, as shown in figure [3j 

Figure [3] shows that the initial amount of entanglement does not uniquely determine 
the remaining entanglement after a propagation distance t, as was found for the two- 
dimensional situation [5J. The Bell states of Set 3 maintain a non-zero entanglement 
for longer than the states of Sets 1 and 2. This can be understood as a consequence of 
the fact that there is a stronger coupling in this case between modes having azimuthal 
indices with different magnitudes than between those with identical magnitudes, which 
in turn follows from how the coupling strength between two OAM modes depends on 
the difference between their OAM-values. 
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Normalized propagation distance Normalized propagation distance 

(a) (b) 

Figure 3: Comparison of the tangle evolution of the Bell states in Set 1, Set 2 and Set 
3 for strong turbulence (a) and for weak turbulence (b). 



5.2. Initially maximally entangled states 

Next we investigate the decay of entanglement of initially maximally entangled qutrit 
states. Depending on their decay curves, one can identify two sets of initial qutrit states, 
given by 

State 1 = [exp (i0i) \j, k) + exp (i0 2 ) -k) + exp (i0 3 ) |0, 0)] (18) 
V3 

State 2 = — [exp (10i) |0, -j) + exp (10 2 ) \-j, 0) + exp (k/> 3 ) \ j,j)] , (19) 

where j, k G {—1, 1}, and <f>i, 02, 03, "02 and -03 are arbitrary phases. Calculating 
the tangle for states with the initial conditions given as State 1 in ([TBI) , we find that 
the phases always appear in the combination = 0i + 02 — 203 in the expression of the 
tangle. The expression for the combined phase that maximizes the tangle at t — to is 
given by 

2 H t (t ) #i(* ) 



= 0o P t (to) = arctan 



(20) 



H^to) 2 - H r (to) 2 . 

while = op t (^o) ~ 7T minimizes the tangle. Note that, although the combined phase 
that minimizes or maximizes the tangle depends on the propagation distance, this 
dependence is very small in the case of strong turbulence, because the entanglement 
decays quickly, for t < 1/3. As a result we can use the small t approximation, which 
implies that 0^0 and ~ n, respectively, maximizes and minimizes the tangle. 
We'll denote the cases with minimum and maximum tangle as State la and State lb, 
respectively. In the case where the initial state is given by State 2 the tangle is totally 
independent of the phases ipi, ip 2 and -03- 

The decay curves of State la, State lb and State 2 are shown in figure HI both 
in weak and strong turbulence. One finds that states represented by State 2 are more 
robust than those represented by State 1. This again follows from the stronger coupling 
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Figure 4: Tangle evolution of the initially maximally entangled state in strong turbulence 
(a) and in weak turbulence (b), both showing three curves: State la, State lb and State 
2. In strong turbulence State la and State lb are obtained for = it and = 0, 
respectively, and in weak turbulence State la and State lb are obtained for = OAw 
and = — 0.67T, respectively. 



between differing magnitudes of azimuthal indices compared to the coupling strength 
between modes with identical magnitudes of azimuthal indices. 

Comparing the curves for strong turbulence (figure Ha| to the ones for weak 
turbulence (figure l4bl) . one can see that the entanglement of states for larger propagation 
distances decays comparatively faster — if we would have used the approximations for 
Z(t), H r (t) and H{(t), given in (J9j), ( TTU1) and ( TTTT) . respectively, to calculate the tangle 
evolution for weak turbulence (i.e. when a = 0.25), we would have found that the 
entanglement was maintained much longer than it does for the exact expressions. This 
can be seen from the fact that the tangle for weak turbulence decays nearly linearly all 
the way to zero entanglement, whereas for strong turbulence it decays linearly up to 
t 0.005, followed by a nonlinear slower decay. Thus, if one can operate in a scenario 
where a and have larger values, then one would have robust states for comparatively 
larger propagation distances. 

5.3. Most robust entangled state 

Upon comparing the curves for initially maximally entangled qutrit states in figure H] 
to those for the Bell states of the two-dimensional subsystems in figure [3J we find that 
the maximally entangled qutrit states are not the most robust qutrit states. The Bell 
states of Set 3 maintain their entanglement longer than any of the initially maximally 
entangled qutrit states in ffl8|) and ffl9l) . even though the initial entanglement of the Bell 
states is smaller. As a result, if one wants to find the most robust entangled state for 
a specific t, it is not enough to consider only those states that are initially maximally 
entangled. 
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Figure 5: Optimization parameters Hi(t) for State 1 and ^(t) for State 2 that gives 
the most robust entangled states \ni{t)) and |/x 2 (t)) as a function of the normalized 
propagation distance t, in strong turbulence (a) and weak turbulence (b). 

To find the most robust states we optimize the tangle at a particular point t = to, for 
the general initial pure state, using the parameterization from ((3]). This optimization 
process is similarly to [8]. From the optimization results, expressed in terms of the 
parameters fc a , k h and g a , we found that all the states have one of the following 

two forms 

lMi> = -j= cos(^i) exp(i0i) \j, k) + -j= cos(//i) exp(i0 2 ) —k) 

+ sin(// 1 )exp(i0 3 )|O,O) (21) 
|/i 2 ) = -j= cos(/i 2 ) exp(i^i) |0, -j) + -j= cos(/i 2 ) exp(k/> 2 ) \-j, 0) 

+ sin(// 2 ) exp(f0 3 ) \j,j) , (22) 

where we introduced the optimization parameters fi\ and /i 2 in addition to the phases 
0i, 02, 03, i>i, ip2 and ip3. These parameters are used to replace those given in ([3]) for 
the sake of simplicity. The optimization parameters /ii and /i 2 determine the degree of 
relative weighting between the maximally entangled states in (Tl8|) and (fl9|) and the Bell 
states of the two-dimensional subspaces given in ( Fl3l) . (fT4l) and (fT5l) . 

The curves of the optimized parameters /ii(t) and /i 2 (t) are shown in figure [SJ for 
both weak and strong turbulence. The fact that the parameters are not independent 
of t means that there does not exist one single state that is the most robust state 
for all propagation distances. Instead there are different states that have the highest 
remaining entanglement at each specific value of t. For t = the optimized states 
|/xi (0)) and |/U 2 (0)) are the maximally entangled states of f|T8|) and ffl9l . respectively. 

Substituting fT20|) and the optimal values for fii and /x 2 at t = 0.03 (t = 1.6) 
for strong (weak) turbulence into ( 12TI) and ( 1221) . we obtain states that maintain a non- 
zero entanglement much longer than the initially maximally entangled states of (fl8l) 
and (|T9l) or the Bell states in ( fl3l) . ( fl4l) and ( fl5l) . The corresponding curves of the 
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Figure 6: Comparison of the tangle evolution of initially maximally entangled states 
(State 1 and State 2) to those of the optimized states (\^i(t)) and \fi 2 (t))) for strong 
turbulence optimized at t = 0.03 (a) and for weak turbulence optimized at t = 1.6 
(b). 

tangle are shown in figure El Note that State 1, which is less robust than State 2, can 
be optimized to produce state that is more robust for large propagation distances 
than the optimized state \^2)- It turns out that, when optimized, is the most 
robust qutrit state within our three-dimensional Hilbert space for propagation through 
turbulence. 

6. Summary and conclusions 

We have shown that for three-dimensional systems, the initial amount of entanglement 
does not uniquely define the amount of entanglement left after a non-zero propagation 
through a turbulent medium. This highlights the fact that the upper bound for the 
entanglement evolution [101 [11] has to be treated with caution. The reason is that in our 
three-dimensional Hilbert space the coupling strength between differing magnitudes of 
azimuthal indices is different from the coupling strength between identical magnitudes. 
This follows from our choice of Hilbert space and the fact that the coupling strength 
decreases as the difference between the azimuthal modal indices increases. 

By optimizing the tangle, starting from the most general initial pure state, we 
found that neither the initially maximally entangled states, nor one of the Bell states 
of the two-dimensional subsystems are the most robust states. Instead it is a coherent 
superposition between both. For each propagation distance one can find a specific initial 
entangled state that retains the most entanglement up to that propagation distance. 
Thus, there is not one single most robust state for all distances. 

The effect of the optimization of the robustness of the qutrit states is still rather 
small (it only increases the propagation distance by a few percent). However, it 
is reasonable to expect that the effect of such an optimization would become more 
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significant for higher dimensional Hilbert spaces, making it worth doing the optimization 
to find robust quantum states in higher dimensional Hilbert spaces and thereby improve 
the performance of quantum communication systems. As a result we believe that 
the optimization of the robustness of higher dimensional quantum states will have a 
significant impact in practical free-space quantum communication systems. 
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